# Volha Charnysh
# Explaining outgroup bias in weak states: Religion and legibility in the 1891-92 Russian famine” (World Politics)
# Replication for panel analysis of mortality and natality during the famine
# R version 3.6.3 
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.6
# ----------

rm(list = ls())
### Install missing packages

ipak <- function(pkg) {
	new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
	if (length(new.pkg)) 
		install.packages(new.pkg, dependencies = TRUE)
	sapply(pkg, require, character.only = TRUE)
}

packages <- c("foreign",  "glmnet", "lfe", "fastDummies", "plyr", "dplyr", "broom", "ggplot2", "RColorBrewer", "reshape2", "rgdal", "geosphere", "RCurl", "texreg", "gplots", "tidyr", "plm", "stargazer", "lmtest", "multiwayvcov", "data.table", "foreign", "reshape", "olsrr")

ipak(pkg = packages) 

#Read in panel dataset: 
dat <- read.dta("time-series-dataset.dta")

############################################################################ 
################## DESCRIPTIVE STATISTICS for Table A.1 #################### 
############################################################################ 

stargazer(dat[dat$aidgub==1 & (!is.na(dat$ostatokppLag)),which(names(dat) %in% c("ostatokpp", "est_cdr", "est_cbr")) ], title="Descriptive statistics", digits=2,summary.stat=c("n", "mean","sd",  "min", "max"),covariate.labels=c( "Deaths per 1000 people", "Births per 1000 people","Harvest, puds per capita"))


###########################################################################################
################## FIGURE 2: PLOT DEATH RATES and HARVEST FAILURE OVER TIME   #############
###########################################################################################

#Above and below average Muslim share for death rates:

mean(dat$sh_muslim_1870[dat$aidgub==1], na.rm=TRUE)

pdf("Figure4a.pdf", width=7, height=5) 
plotmeans(est_cdr ~ year, ylab = "Death rates (per 1,000 people)", xlab="", data = subset(dat,sh_muslim_1870>=mean(dat$sh_muslim_1870[dat$aidgub==1], na.rm=TRUE) & aidgub==1), n.label = FALSE, pch = 19, ylim = c(33, 56), xaxs="i")
plotmeans(est_cdr ~ year, data = subset(dat, sh_muslim_1870<mean(dat$sh_muslim_1870[dat$aidgub==1], na.rm=TRUE) & aidgub ==1), n.label = FALSE, add = TRUE, lty = 2, pch = 21, xaxs="i")
legend("topleft", title="Share Muslim at the district level", legend = c("above average (6.4%)", "below average (6.4%)"), 
	lty = 1:2, pch = c(19, 21), cex = 0.8, box.lty = 0)
dev.off() 

#Above and below average Muslim share for birth rates: 

pdf("Figure4b.pdf", width=7, height=5) 
plotmeans(est_cbr ~ year, ylab = "Birth rates (per 1,000 people)",xlab="", data = subset(dat, sh_muslim_1870>=mean(dat$sh_muslim_1870[dat$aidgub==1], na.rm=TRUE) & aidgub ==1), n.label = FALSE, pch = 19, ylim = c(43, 60), xaxs="i")
plotmeans(est_cbr ~ year, data = subset(dat, sh_muslim_1870<mean(dat$sh_muslim_1870[dat$aidgub==1], na.rm=TRUE) & aidgub ==1),  n.label = FALSE, add = TRUE, lty = 2, pch = 21, xaxs="i")
legend("bottomleft", title="Share Muslim at the district level", legend = c("above average (6.4%)", "below average (6.4%)"), 
	lty = 1:2, pch = c(19, 21), cex = 0.8, box.lty = 0)
dev.off() 


######################################################################
######################### Create Functions ###########################
######################################################################


formulas <- function(a, b, c) {
  formulas <- list()
  n <- 0
  for(k in c){
    for (i in a) {
      for (j in b) {
        n <- n + 1
        formulas[n] <- paste(i, j, k)
      }
    }
  }
  return(formulas)
}


run_models <- function(formlist) {
  models <- list()
  n <- 0
  for (i in 1:length(formlist)){
    n <- n + 1
    models[[n]] <- eval(parse(text=formlist[[n]])) 
  }
  return(models)
}

#function to extract coefficients + calculate errors with VCOV
run_coef <- function(output, datalist) {
  stdev <- list() 
  n <- 0
  for (i in 1:length(output)){
    n <- n + 1
    stdev[[n]] <- coeftest(output[[n]], vcov = cluster.vcov(output[[n]], datalist))
  }
  return(stdev)    
}


##################################################################################
################## Table 1: Two-way fixed-effects regression  ####################
##################################################################################

#separate out DVs, rel variables, and controls
dv <- list("lm(est_cdr ~ factor(ter_id) + factor(year) +", "lm(est_cbr ~ factor(ter_id) + factor(year) +")
religion <- list("I(famine*sh_NOTorthodox_1870)", "I(famine*sh_muslim_1870) + I(famine*sh_NOTorthNotMuslim1870)", "I(famine*ShareTurk) + I(famine*ShareNonRusNonTurk)")
controls <- list("+ ostatokppLag + ostatokppLag*famine + UrbanShare1883*famine, data = subset(dat,  dat$aidgub==1))")

#create and run formulas + get coefficients
reg.formula <- formulas(dv, religion, controls)
mod <- run_models(reg.formula)
mod.cl <- run_coef(mod, dat$ter_id[dat$aidgub==1])

#CREATING TABLE 1: Religion, language, and district-level mortality and natality during the famine. 
tableMain <- texreg(
  l = list(mod[[1]], mod[[2]], mod[[3]], mod[[4]], mod[[5]], mod[[6]]),
  override.se = list(mod.cl[[1]][,2], mod.cl[[2]][,2], mod.cl[[3]][,2], mod.cl[[4]][,2], mod.cl[[5]][,2], mod.cl[[6]][,2]),
  override.pval = list(mod.cl[[1]][,4], mod.cl[[2]][,4], mod.cl[[3]][,4], mod.cl[[4]][,4], mod.cl[[5]][,4], mod.cl[[6]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)", "(5)", "(6)"),
  custom.coef.names=c("Famine X non-Orthodox", "Harvest per capita (lag)", "Famine X Harvest per capita (lag)", "Famine X Share Urban", "Famine X Muslims","Famine X Other Non-Orthodox", "Famine X Turkic","Famine X Other Non-Russian"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = "(Intercept)|ter_id|year",
  reorder.coef = c(1,5:8,2:4),
  include.rsquared = FALSE, 
  caption="Religion, language, and district-level mortality and natality during the famine.",
  caption.above = TRUE,
)


##################################################################################
############## Table A3 (Appendix) with famine coded 1892 AND 1893  ##############
##################################################################################

#separate out DVs, rel variables, and controls
dv.1a <- list("lm(est_cdr ~ factor(ter_id) + factor(year) +", "lm(est_cbr ~ factor(ter_id) + factor(year) +")
religion.1a <- list("I(famine2*sh_NOTorthodox_1870)", "I(famine2*sh_muslim_1870) + I(famine2*sh_NOTorthNotMuslim1870)", "I(famine2*ShareTurk) + I(famine2*ShareNonRusNonTurk)")
controls.1a <- list("+ ostatokppLag + ostatokppLag*famine2 + UrbanShare1883*famine2, data = subset(dat,  dat$aidgub==1))")

#create and run formulas + get coefficients
reg.formula.1a <- formulas(dv.1a, religion.1a, controls.1a)
mod.1a <- run_models(reg.formula.1a)
mod.1a.cl <- run_coef(mod.1a, dat$ter_id[dat$aidgub==1])

#CREATING TABLE A.3 in the Appendix.
tableMain.3a <- texreg(
  l = list(mod.1a[[1]], mod.1a[[2]], mod.1a[[3]], mod.1a[[4]], mod.1a[[5]], mod.1a[[6]]),
  override.se = list(mod.1a.cl[[1]][,2], mod.1a.cl[[2]][,2], mod.1a.cl[[3]][,2], mod.1a.cl[[4]][,2], mod.1a.cl[[5]][,2], mod.1a.cl[[6]][,2]),
  override.pval = list(mod.1a.cl[[1]][,4], mod.1a.cl[[2]][,4], mod.1a.cl[[3]][,4], mod.1a.cl[[4]][,4], mod.1a.cl[[5]][,4], mod.1a.cl[[6]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)", "(5)", "(6)"),
  custom.coef.names=c("Famine X non-Orthodox", "Harvest per capita (lag)", "Famine X Harvest per capita (lag)", "Famine X Share Urban", "Famine X Muslims","Famine X Other Non-Orthodox", "Famine X Turkic","Famine X Other Non-Russian"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = "(Intercept)|ter_id|year",
  reorder.coef = c(1,5:8,2:4),
  include.rsquared = FALSE, 
  caption="Religion, language, and district-level mortality and natality during the famine.The famine dummy was coded one for both 1892 and 1893.",
  caption.above = TRUE,
)

############################################################################################################################
################ Table A.4 with Conley standard errors was produced in Stata. See code in Conley_SE.do  ####################
############################################################################################################################


############################################################################################################################
################ Table A.5: Two-way fixed-effects regression with additional controls suggested by reviewer ################
############################################################################################################################

dat$AvgRecrHeightM <- dat$AvgRecrHeight/100 #change recruit heigh from cm to meters.
  
dv.A5 <- list("lm(est_cdr ~ factor(ter_id) + factor(year) +", "lm(est_cbr ~ factor(ter_id) + factor(year) +")
religion.A5 <- list("I(famine*sh_NOTorthodox_1870)", "I(famine*sh_muslim_1870) + I(famine*sh_NOTorthNotMuslim1870)", "I(famine*ShareTurk) + I(famine*ShareNonRusNonTurk)")
controls.A5 <- "+ ostatokppLag + ostatokppLag*famine + UrbanShare1883*famine + lnrailDist*famine + serfperc1*famine + AvgRecrHeightM*famine + DvorsWithHorsesPct*famine + nobilityper1000pop77*famine + captPerArea*famine + lndistStPetKm*famine + PopDensity*famine, data = subset(dat,  dat$aidgub==1))"

#create and run formulas + get coefficients
reg.formula.A5 <- formulas(dv.A5, religion.A5, controls.A5)
mod.A5 <- run_models(reg.formula.A5)
mod.cl.A5 <- run_coef(mod.A5, dat$ter_id[dat$aidgub==1])

#CREATING TABLE: 
tableA5 <- texreg(
  l = list(mod.A5[[1]], mod.A5[[2]], mod.A5[[3]], mod.A5[[4]], mod.A5[[5]], mod.A5[[6]]),
  override.se = list(mod.cl.A5[[1]][,2], mod.cl.A5[[2]][,2], mod.cl.A5[[3]][,2], mod.cl.A5[[4]][,2], mod.cl.A5[[5]][,2], mod.cl.A5[[6]][,2]),
  override.pval = list(mod.cl.A5[[1]][,4], mod.cl.A5[[2]][,4], mod.cl.A5[[3]][,4], mod.cl.A5[[4]][,4], mod.cl.A5[[5]][,4], mod.cl.A5[[6]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)", "(5)", "(6)"),
  custom.coef.names=c("Famine X non-Orthodox", "Harvest per capita (lag)", "Famine X Harvest per capita (lag)", "Famine X Share Urban", "Famine X Ln(Distance to railway)", "Famine X Serfdom (share)", "Famine X Avg. Recruit Heights", "Famine X Horses per household", "Famine X Noble Landowners" ,"Famine X Land Captains", "Famine X Ln(Distance to St. Petersburg)", "Famine X Population Density", "Famine X Muslims","Famine X Other Non-Orthodox", "Famine X Turkic", "Famine X Other Non-Russian"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = "(Intercept)|ter_id|year",
  reorder.coef = c(1,13:16, 2:12),
  include.rsquared = FALSE
)


############################################################################################################
################## Table A.6: Two-way fixed-effects regression with geographic controls ####################
############################################################################################################

#change controls and religion variables
dv.A6 <- list("lm(est_cdr ~ factor(ter_id) + factor(year) +", "lm(est_cbr ~ factor(ter_id) + factor(year) +")
religion.A6 <- list("I(famine*sh_muslim_1870) + I(famine*sh_NOTorthNotMuslim1870)")
controls.A6 <- list("+ ostatokppLag + ostatokppLag*famine + UrbanShare1883*famine +I(lat*famine), data = subset(dat,  dat$aidgub==1))", "+ ostatokppLag + ostatokppLag*famine + UrbanShare1883*famine + I(long*famine), data = subset(dat,  dat$aidgub==1))", "+ ostatokppLag + ostatokppLag*famine + UrbanShare1883*famine + I(long*lat*famine), data = subset(dat,  dat$aidgub==1))")

#create and run formulas + get coefficients
reg.formula.A6 <- formulas(dv.A6, religion.A6, controls.A6)
mod.A6 <- run_models(reg.formula.A6)
mod.cl.A6 <- run_coef(mod.A6, dat$ter_id[dat$aidgub==1])

#CREATING TABLE: 
tableGeoCovs <- texreg(
  l = list(mod.A6[[1]], mod.A6[[3]], mod.A6[[5]], mod.A6[[2]], mod.A6[[4]], mod.A6[[6]]),
  override.se = list(mod.cl.A6[[1]][,2], mod.cl.A6[[3]][,2], mod.cl.A6[[5]][,2], mod.cl.A6[[2]][,2], mod.cl.A6[[4]][,2], mod.cl.A6[[6]][,2]),
  override.pval = list(mod.cl.A6[[1]][,4], mod.cl.A6[[3]][,4], mod.cl.A6[[5]][,4], mod.cl.A6[[2]][,4], mod.cl.A6[[4]][,4], mod.cl.A6[[6]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)", "(5)", "(6)"),
  custom.coef.names=c("Famine X Muslims","Famine X Other Non-Orthodox", "Harvest per capita (lag)", "Famine X Latitude", "Famine X Harvest per capita (lag)", "Famine X Share Urban", "Famine X Longitude", "Latitude X Longitude X Famine"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = "(Intercept)|ter_id|year",
  include.rsquared = FALSE
)

########################################################################################################
############## Table A.7: ESTIMATE COEFFICIENT ON "sh_muslim_1870" separately for each year ############
########################################################################################################

dv.A7 <- list("lm(est_cdr ~  factor(ter_id) + factor(year) +", "lm(est_cbr ~  factor(ter_id) + factor(year) +")
religion.A7 <- list("factor(year)*sh_muslim_1870 + famine*sh_NOTorthNotMuslim1870", "factor(year)*ShareTurk + famine*ShareNonRusNonTurk")
controls.A7 <- list("+ ostatokppLag + I(ostatokppLag*famine) + I(UrbanShare1883*famine), subset(dat,  dat$aidgub==1))")

#create and run formulas + get coefficients
reg.formula.A7 <- formulas(dv.A7, religion.A7, controls.A7)
mod.A7 <- run_models(reg.formula.A7)
mod.cl.A7 <- run_coef(mod.A7, dat$ter_id[dat$aidgub==1])

#CREATING TABLE A.7: Note the table was formatted manually with some dummies removed to fit on the page in the Appendix.
tableMainA7 <- texreg(
  l = list(mod.A7[[1]], mod.A7[[2]], mod.A7[[3]], mod.A7[[4]]),
  override.se = list(mod.cl.A7[[1]][,2], mod.cl.A7[[2]][,2], mod.cl.A7[[3]][,2], mod.cl.A7[[4]][,2]),
  override.pval = list(mod.cl.A7[[1]][,4], mod.cl.A7[[2]][,4], mod.cl.A7[[3]][,4], mod.cl.A7[[4]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)"),
  custom.coef.names=c("year 1890", "year 1891", "year 1892", "year 1893", "year 1894", "year 1895", "year 1896", "Harvest per capita (lag)","Famine X Harvest per capita (lag)", "Famine X Share Urban", 
  "year 1890 X Muslims", "year 1891 X Muslims","year 1892 X Muslims", "year 1893 X Muslims", "year 1894 X Muslims", "year 1895 X Muslims", "year 1896 X Muslims", "Famine X Other non-Orthodox", "year 1890 X Share Turkic", "year 1891 X Share Turkic", "year 1892 X Share Turkic", "year 1893 X Share Turkic", "year 1894 X Share Turkic", "year 1895 X Share Turkic", "year 1896 X Share Turkic", "Famine X Other non-Russian"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = c("(Intercept)|ter_id"),
  include.rsquared = FALSE, 
  caption="Religion, language, and district-level mortality and natality during the famine. Famine dummy is equal to one for 1892",
  caption.above = TRUE,
)

##################################################################################### 
################## Using regularization to pick relevant controls  ##################
#####################################################################################

#First create a subset of the data that excludes NAs 
dat3 <- as.data.frame(dat)
dat3 <- subset(dat3, year>1888  & aidgub==1)

vars<-c("est_cdr","est_cbr", "famine", "sh_NOTorthodox_1870", "sh_muslim_1870", "sh_NOTorthNotMuslim1870", "ShareTurk", "ShareNonRusNonTurk", "ter_id", "year", "ostatokppLag", "UrbanShare1883", "lnrailDist", "serfperc1", "DvorsWithHorsesPct", "nobilityper1000pop77","AvgRecrHeightM", "captPerArea", "lndistStPetKm","PopDensity")

dat3 <- dat3[,vars]
dat3 <- na.omit(dat3) 

#create dummies from X for estimation
dat3 <- fastDummies::dummy_cols(dat3, select_columns = c("ter_id", "year")) 

#cluster by district
cluster <- factor(dat3$ter_id)

#There are two main outcomes, estimated death rate and estimated birth rate; create list
cdr<-dat3$est_cdr
cbr<-dat3$est_cbr

outcomes <- list(cdr, cbr)

#There are different treatment measures, share non-Orthodox x famine dummy, Muslim x famine dummy, and Turk x famine dummy; create list.
treat1<-data.matrix(dat3$sh_NOTorthodox_1870*dat3$famine)
treat2<-data.matrix(dat3$sh_muslim_1870*dat3$famine)
treat3<-data.matrix(dat3$ShareTurk*dat3$famine)

treatment <- list(treat1, treat2, treat3)

#create data matrix 
ds.matrix1 <- bind_cols(dat3$ostatokppLag, dat3$ostatokppLag*dat3$famine, dat3$UrbanShare1883*dat3$famine, dat3$lnrailDist*dat3$famine, dat3$serfperc1*dat3$famine, dat3$AvgRecrHeightM*dat3$famine, dat3$DvorsWithHorsesPct*dat3$famine, dat3$nobilityper1000pop77*dat3$famine, dat3$captPerArea*dat3$famine, dat3$lndistStPetKm*dat3$famine, dat3$PopDensity*dat3$famine, dat3[,21:244]) 

ds.matrix2 <- bind_cols(dat3$sh_NOTorthNotMuslim1870*dat3$famine, dat3$ostatokppLag, dat3$ostatokppLag*dat3$famine, dat3$UrbanShare1883*dat3$famine, dat3$lnrailDist*dat3$famine, dat3$serfperc1*dat3$famine, dat3$AvgRecrHeightM*dat3$famine, dat3$DvorsWithHorsesPct*dat3$famine, dat3$nobilityper1000pop77*dat3$famine, dat3$captPerArea*dat3$famine, dat3$lndistStPetKm*dat3$famine, dat3$PopDensity*dat3$famine, dat3[,21:244]) 

ds.matrix3 <- bind_cols(dat3$ShareNonRusNonTurk*dat3$famine, dat3$ostatokppLag, dat3$ostatokppLag*dat3$famine, dat3$UrbanShare1883*dat3$famine, dat3$lnrailDist*dat3$famine, dat3$serfperc1*dat3$famine, dat3$AvgRecrHeightM*dat3$famine, dat3$DvorsWithHorsesPct*dat3$famine, dat3$nobilityper1000pop77*dat3$famine, dat3$captPerArea*dat3$famine, dat3$lndistStPetKm*dat3$famine, dat3$PopDensity*dat3$famine, dat3[,21:244]) 

colnames(ds.matrix1)<-c("ostatokppLag", "ostatokppLagXfamine", "UrbanShare1883Xfamine", "lnrailDistXfamine", "serfperc1Xfamine","AvgRecrHeightMXfamine","DvorsWithHorsesPctXfamine", "nobilityper1000pop77Xfamine", "captPerAreaXfamine", "lndistStPetKmXfamine", "PopDensityXfamine", names(dat3)[21:244])

colnames(ds.matrix2)<-c("sh_NOTorthNotMuslim1870Xfamine", "ostatokppLag", "ostatokppLagXfamine", "UrbanShare1883Xfamine", "lnrailDistXfamine", "serfperc1Xfamine","AvgRecrHeightMXfamine", "DvorsWithHorsesPctXfamine", "nobilityper1000pop77Xfamine", "captPerAreaXfamine", "lndistStPetKmXfamine", "PopDensityXfamine",  names(dat3)[21:244])

colnames(ds.matrix3)<-c("ShareNonRusNonTurkXfamine" , "ostatokppLag", "ostatokppLagXfamine", "UrbanShare1883Xfamine", "lnrailDistXfamine", "serfperc1Xfamine", "AvgRecrHeightMXfamine","DvorsWithHorsesPctXfamine", "nobilityper1000pop77Xfamine", "captPerAreaXfamine", "lndistStPetKmXfamine", "PopDensityXfamine", names(dat3)[21:244])

ds.matrix1 <- data.matrix(ds.matrix1) 
ds.matrix2 <- data.matrix(ds.matrix2) 
ds.matrix3 <- data.matrix(ds.matrix3) 

ds.matrices <- list(ds.matrix1, ds.matrix2, ds.matrix3)


#################################################################
########### Table A.8: USING DOUBLE-SELECTION METHOD ############
#################################################################

#Double-selection models outcome and treatment separately and uses only variables selected by BOTH models

double_selection_fn <- function(treatment, outcome, datamatrix) {
  models <- list()
  n <- 0
  set.seed(1000)
  
  for (i in 1:length(outcome)){
    for (j in 1:length(treatment)) {
      
      nFolds <- 10
      foldid <- sample(rep(seq(nFolds), length.out = nrow(treatment[[i]])))
      
      n <- n + 1
      
      a <-  cv.glmnet(x = cbind(treatment[[j]], datamatrix[[j]]), y = outcome[[i]], nfolds = nFolds, foldid = foldid)
      b <- as.matrix(coef(a)[-c(1:2)]) != 0
      
      c <- cv.glmnet(x = datamatrix[[j]], y = treatment[[j]],  nfolds = nFolds, foldid = foldid)
      d <- as.matrix(coef(c)[-1]) != 0
      
      models[[n]] <- lm(outcome[[i]] ~ treatment[[j]] + datamatrix[[j]][, b|d])
    }
  }
  return(models)
}

double_selection_fn_SE <- function(cvglm, cluster) {
  SE_models <- list()
  n <- 0
  for (i in 1:length(cvglm)){
    n <- n + 1
    
    SE_models[[n]] <- coeftest(cvglm[[i]], vcov = cluster.vcov(cvglm[[i]], cluster))
  }
  return(SE_models)
}

#store double selection in 'ds' and errors in 'ds.SE'
ds <- double_selection_fn(treatment, outcomes, ds.matrices)
ds.SE <- double_selection_fn_SE(ds, cluster)

#create table with results
#Note treatment here varies by the model: Share non-Orthodox (1, 4), share Muslim (2, 5), Share Turkic (3,6) 
tableDS <- texreg(
  l = list(ds[[1]], ds[[2]], ds[[3]], ds[[4]], ds[[5]], ds[[6]]),
  override.se = list(ds.SE[[1]][,2], ds.SE[[2]][,2], ds.SE[[3]][,2], ds.SE[[4]][,2], ds.SE[[5]][,2], ds.SE[[6]][,2]),
  override.pval = list(ds.SE[[1]][,4], ds.SE[[2]][,4], ds.SE[[3]][,4], ds.SE[[4]][,4], ds.SE[[5]][,4], ds.SE[[6]][,4]),
  custom.model.names = c("(1)","(2)","(3)","(4)", "(5)", "(6)"),
  custom.coef.names=c("Treatment *", "Harvest per capita (lag)", "Famine X Harvest per capita (lag)", "Famine X Share Urban", "Famine X Ln(Distance to railway)", "Famine X Serfdom (share)","Famine X Avg. Recruit Heights", "Famine X Horses per household", "Famine X Noble Landowners", "Famine X Land Captains", "Famine X Ln(Distance to St. Petersburg)", "Famine X Population Density",  "Famine X Other non-Orthodox", "Famine X Other non-Russian"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:mortalityTSU",
  omit.coef = "(Intercept)|ter_id|year",
  include.rsquared = FALSE, 
  caption="Table A.8: Religion, language, and district-level mortality and natality during the famine.Famine dummy is coded one for 1892. Estimates are based on the double-selection methodfor including covariates proposed by Chernozhukov et al. (2018).",
  caption.above=TRUE,
)